
#---- Main analysis ------#
# This file conducts the descriptive and multivariate analyses as well as the appendix as shown in
# "Supply-side dynamics of group appeals: How dominance affects 
# parties' choice between symbolic and policy-based appeals" (2025)

# ---- 0. Set working directory (include correct path, create a "Data" and "Figures" folder) ----
setwd("")

# ---- 1. Load packages needed -----
library(tidyverse)
library(ggpubr)
library(openxlsx)
library(gridExtra)
library(stargazer)

#
# ---- 2. Load data ----
#


# ---- 2.1 Load analysis dataset -----
groupappeals_fullset_v1 <- readRDS("Data/groupappeals_fullset_v1.rds") 

# ----- 2.2 Load dominance data -----
load("Data/dominance_since1990.RData") 


##
#---- 3. Descriptive Analysis------
##


#
# ----- 3.1 Groups by frequency and % policy-based (bubble plot) (Figure 2) ------
#

# share of group appeals that are accompanied by policy appeals

# prepare data
groupappeals_fullset_v1$subst_policy <- as.factor(groupappeals_fullset_v1$subst_policy)

# create new df
groups <- groupappeals_fullset_v1

# change count variable to one for each group (because the number of mentions 
# per tweet is not relevant)
for (i in 11:35) {
  
  groups[i] <- ifelse(groups[i] > 0,1,0)
}

# aggregate group references
groups <- groups %>% 
  group_by(screen_name, partei, subst_policy) %>%
  summarize(across(students:landlords, sum)) 

# further adjustments
groups$academics <- NULL
groups$screen_name <- NULL
groups$partei <- as.character(groups$partei)

# calculate group totals (##NEW##)
groups_rel <- groups %>%
  dplyr::group_by(partei, subst_policy) %>%
  summarize(across(students:landlords,sum, na.rm = TRUE))

# absolute values and totals
groups_rel <- groups_rel %>%
  group_by(partei) %>%
  mutate(total = sum(across(students:landlords)))

# relative
groups_rel  <- groups_rel  %>%
  group_by(partei, subst_policy) %>%
  summarize(round(across(students:landlords, ~ . /total), digits = 3))


# Helper function for the plot
alpha_values <- function(vec){
  
  vec[vec == 0] <- -1
  vec[vec > 0] <- 0
  vec[vec == -1] <- 100
  return(vec)
  
}


# Get the column names in descending alphabetical order
column_names <- names(groups_rel[3:26])
sorted_column_names <- sort(column_names, decreasing = TRUE)

# Create a new dataframe with the desired column order
groups_rel <- cbind(groups_rel[c(1:2)], groups_rel[sorted_column_names])

# change values where no group appeal was issues to zero
groups_rel[3:26][groups_rel[3:26] == 0] <- NA

# exclude free voters
groups_rel <- groups_rel %>% filter(partei != "FW")

# create ggplot object
results <- ggplot()

# append values for each group and mode (symbolic/policy-based) separately
results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(3-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,3])),
                                    color = "black", pch = 21
)


results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(3-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,3])),
                                color = "black", pch = 21
)



results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(4-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,4])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(4-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,4])),
                                color = "black", pch = 21
)



results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(5-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,5])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(5-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,5])),
                                color = "black", pch = 21
)



results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(6-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,6])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(6-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,6])),
                                color = "black", pch = 21
)




results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(7-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,7])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(7-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,7])),
                                color = "black", pch = 21
)


results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(8-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,8])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(8-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,8])),
                                color = "black", pch = 21
)



results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(9-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,9])),
                                color = "black", pch = 21
)


results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(9-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,9])),
                                color = "black", pch = 21
                                
)


results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(10-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,10])),
                                color = "black", pch = 21
)

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(10-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,10])),
                                color = "black", pch = 21 
)


results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(11-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,11])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(11-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,11])),
                                color = "black", pch = 21 
)


results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(12-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,12])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(12-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,12])),
                                color = "black", pch = 21
)



results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(13-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,13])),
                                color = "black", pch = 21 
)


results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(13-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,13])),
                                color = "black", pch = 21 
)


results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(14-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,14])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(14-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,14])),
                                color = "black", pch = 21 
)



results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(15-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,15])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(15-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,15])),
                                color = "black", pch = 21 
)


results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(16-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,16])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(16-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,16])),
                                color = "black", pch = 21 
)


results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(17-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,17])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(17-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,17])),
                                color = "black", pch = 21 
)


results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(18-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,18])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(18-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,18])),
                                color = "black", pch = 21 
)


results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(19-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,19])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(19-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,19])),
                                color = "black", pch = 21 
)


results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(20-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,20])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(20-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,20])),
                                color = "black", pch = 21 
)



results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(21-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,21])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(21-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,21])),
                                color = "black", pch = 21 
)



results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(22-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,22])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(22-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,22])),
                                color = "black", pch = 21 
)


results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(23-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,23])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(23-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,23])),
                                color = "black", pch = 21 
)

results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(24-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,24])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(24-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,24])),
                                color = "black", pch = 21 
)


results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(25-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,25])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(25-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,25])),
                                color = "black", pch = 21 
)


results <- results + geom_point(aes(x = c(1:6),
                                    y = rep(26-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 0]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 0,][,26])),
                                color = "black", pch = 21
) 

results <- results + geom_point(aes(x = c(10:15),
                                    y = rep(26-2, 6),
                                    fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                                    size = unlist(groups_rel[groups_rel$subst_policy == 1,][,26])),
                                color = "black", pch = 21 
)


# create party labels
party_labels <- data.frame(
  x = c(1:6, 10:15),
  party = rep(unique(groups_rel$partei), 2)
)


# Define custom party colors
party_colors <- c(
  "CDU" = "black",
  "SPD" = "darkred",
  "Grüne" = "darkgreen",
  "FDP" = "yellow",
  "AfD" = "blue",
  "Linke" = "pink"
)

# create final plot
bubble_plot <- results +
  geom_point(aes(x = c(10:15),
                 y = rep(23-2, 6),
                 fill = unlist(groups_rel$partei[groups_rel$subst_policy == 1]),
                 size = unlist(groups_rel[groups_rel$subst_policy == 1,][,23])),
             color = "black", pch = 21) +
  scale_fill_manual(
    name = "Party",
    values = party_colors
  ) +
  theme_test() +
  theme(axis.title = element_blank()) +
  theme(axis.title = element_blank(),
        text = element_text(family = "Arial", size = 14),
        legend.position = "none")

# Add party labels as a secondary x-axis
bubble_plot +
  geom_text(data = party_labels, aes(x = x, y = 25, label = party), 
            vjust = -0.5, hjust = 0.5, size = 4,
            color = "black", fontface= "bold") +
  scale_y_continuous(breaks = c(1:24), labels = colnames(groups_rel)[3:26]) +
  scale_x_continuous(breaks = c(4, 12), labels = c("Symbolic", "Policy-based"), 
                     name = "") +
  theme(
    panel.grid.major.y = element_line(color = "gray", size = 0.3),  
    panel.grid.minor.y = element_blank()  
  )


#
# 3.2 Policy-based appeals (subst. = substantial) per party (Figure 3)
#

#prepare dataset (subst_perparty) with shares of policy-based appeals 
groupappeals_fullset_v1$subst_policy <- as.numeric(groupappeals_fullset_v1$subst_policy)
groupappeals_fullset_v1$subst_policy[groupappeals_fullset_v1$subst_policy == 1] <- 0
groupappeals_fullset_v1$subst_policy[groupappeals_fullset_v1$subst_policy == 2] <- 1

subst_perparty <- groupappeals_fullset_v1 %>%
  group_by(partei,land) %>%
  mutate(sum_appeal = n()) %>%
  mutate(sum_subst = sum(subst_policy)) %>%
  mutate(prop_subst = sum(subst_policy)/sum_appeal) %>%
  mutate(landesverband = str_c(partei,"_",land))

subst_perparty <- subst_perparty %>%
  dplyr::select(landesverband,sum_appeal, sum_subst, prop_subst) %>%
  distinct()

# merge share of policy-based appeals to dominance data 
subst_perparty <- merge(subst_perparty, dominance_since1990, by = c("partei", 
                                                                    "land"))

# convert dominance to years
subst_perparty$dominance_year <- subst_perparty$dominance_total / 365

# Create plot in figure 3
plot_substperparty <- ggplot(subst_perparty, aes(x = dominance_year, y = prop_subst)) +
  geom_point(aes(color = partei), size = 2, alpha = 0.7) +  # Color dots by party
  geom_smooth(method = lm, se = TRUE, color = "black") +  # Single regression line
  xlab("Dominance in years") + 
  ylab("Share of policy-based group appeals") +
  theme_bw() +
  theme(axis.text = element_text(family = "Arial", size = 14),
        axis.title = element_text(family = "Arial", size = 16),
        legend.position = "bottom") +
  scale_color_manual(values = c("CDU" = "black",
                                "SPD" = "darkred",
                                "Grüne" = "darkgreen",
                                "FDP" = "yellow",
                                "AfD" = "blue",
                                "Linke" = "pink"
  )) +
  labs(color = "Party")

# save to reimport in sentiment script, combining plots
# save(subst_perparty, file = "Data/subst_perparty.RData")





#
#
# ----- tabulated data ------
#

# Create the table
table_data <- table(groupappeals_fullset_v1$dominance_bin, groupappeals_fullset_v1$subst_policy)

# Convert to percentages by row
percentages <- prop.table(table_data, margin = 1) * 100

# Print the result
percentages


#
#---- 4. Multivariate Analysis------
#

#
# ---- 4.1. Dataset preparations for the regression analyses -----
#

# convert variables to appropriate classes
groupappeals_fullset_v1$subst_policy <- as.factor(groupappeals_fullset_v1$subst_policy)
groupappeals_fullset_v1$dominance_bin <- as.factor(groupappeals_fullset_v1$dominance_bin)
groupappeals_fullset_v1$partei <- as.factor(groupappeals_fullset_v1$partei)
groupappeals_fullset_v1$gov <- as.factor(groupappeals_fullset_v1$gov)
groupappeals_fullset_v1$MP_atm <- as.factor(groupappeals_fullset_v1$MP_atm)

# convert electoral proximity variables from days to months
groupappeals_fullset_v1$proximity_nextelection <- groupappeals_fullset_v1$proximity_nextelection / 31
groupappeals_fullset_v1$proximity_btw <- groupappeals_fullset_v1$proximity_btw / 31

# dominance yearly (for all dominance variables)
groupappeals_fullset_v1$dominance_year <- 
  groupappeals_fullset_v1$dominance / 365

groupappeals_fullset_v1$dominance_total_year <- 
  groupappeals_fullset_v1$dominance_total / 365

groupappeals_fullset_v1$MP_year <- 
  groupappeals_fullset_v1$MP / 365

groupappeals_fullset_v1$MP_total_year <- 
  groupappeals_fullset_v1$MP_total / 365

groupappeals_fullset_v1$dominance_before2000_year <- 
  groupappeals_fullset_v1$dominance_before2000 / 365

groupappeals_fullset_v1$dominance_00to10_year <- 
  groupappeals_fullset_v1$dominance_00to10 / 365

groupappeals_fullset_v1$dominance_after2010_year <- 
  groupappeals_fullset_v1$dominance_after2010 / 365

groupappeals_fullset_v1$dominance_wgt_year <- 
  groupappeals_fullset_v1$dominance_wgt / 365

groupappeals_fullset_v1$MP_wgt_year <- 
  groupappeals_fullset_v1$MP_wgt / 365

# create variable indicating whether a party has ever been the PM party
groupappeals_fullset_v1$MP_bin <- 
  ifelse(groupappeals_fullset_v1$MP_total == 0, 0,1)

# specify reference category for party variable
groupappeals_fullset_v1$partei <- as.factor(groupappeals_fullset_v1$partei)
groupappeals_fullset_v1$partei <- relevel(
  groupappeals_fullset_v1$partei, 
  ref = 2)


# measure proximity instead of distance
groupappeals_fullset_v1$proximity_nextelection <- groupappeals_fullset_v1$proximity_nextelection*(-1)
groupappeals_fullset_v1$proximity_btw <- groupappeals_fullset_v1$proximity_btw*(-1)

# create year variable
groupappeals_fullset_v1$year <- substr(groupappeals_fullset_v1$month, 1, 4)

#
# 4.2 Calculate logit models
#

logit_m1 <- glm(subst_policy ~ dominance_bin +
                  partei,
                  #as.factor(land),
                data = groupappeals_fullset_v1 , 
                family= "binomial")


logit_m2 <- glm(subst_policy ~ dominance_year +
                  partei,
                  #as.factor(land),
                data = groupappeals_fullset_v1 , 
                family= "binomial")


logit_m3 <- glm(subst_policy ~ dominance_bin*proximity_nextelection +
                  proximity_btw+
                  partei,
                  #gov +
                  #as.factor(land),
                data = groupappeals_fullset_v1 , 
                family= "binomial")

logit_m4 <- glm(subst_policy ~ dominance_year*MP_year + 
                   partei + 
                   gov,
                   #as.factor(land),
                 data = groupappeals_fullset_v1, 
                 family= "binomial")


## Function to compute odds ratios
extract_or <- function(model) {
  exp(coef(model))  # Convert log-odds to odds ratios
}

# Function to compute confidence intervals in OR scale
compute_ci <- function(model) {
  exp(confint(model))  # Convert log-odds CI to odds ratios CI
}

# Function to extract correct p-values
extract_p_values <- function(model) {
  coefs <- summary(model)$coefficients
  coefs[, 4]  # Extracts p-values from the log-odds model
}

# List of models
models <- list(logit_m4, logit_m6, logit_m8, logit_m11)

# Compute odds ratios, confidence intervals, and p-values
or_list <- lapply(models, extract_or)
ci_list <- lapply(models, compute_ci)
p_values_list <- lapply(models, extract_p_values)

# Export models
stargazer(models, 
          type = "html",
          coef = or_list,
          ci = TRUE,
          ci.custom = ci_list,
          p = p_values_list,  
          omit = "partei",  
          digits = 3,
          title = "Odds Ratios with Correct P-Values and Confidence Intervals",
          dep.var.labels.include = FALSE,
          add.lines = list(c("Party dummies", "YES", "YES", "YES", "YES")),
          out = "Figures/regtable_main_mar2025.html")






#
# ---- 4.2 Average marginal effect of dominance (in years) moderated by MP years (Figure 5) -----
#

# calculate marginal effects of the interaction term
margins_interaction <- margins(logit_m11,
                               at = list(MP_year = c(1:30)),
                               variables = c("dominance_year")
)
sum_interaction <- summary(margins_interaction)
sum_interaction$y <- 1:nrow(sum_interaction)

# Plot marginal effects of dominance by PM in years

plot_amesbyMPyears <- ggplot(data = sum_interaction,
                             aes(y = AME, x = y, ymin = lower, ymax = upper,
                                 xmin= 0, xmax = 30)) +
  geom_point(color = "black") +
  geom_errorbar(width = 0) +
  scale_x_continuous(name = "PM in years",
                     breaks=c(0, 5,10,15,20,25,30)
  ) +
  geom_hline(yintercept= 0, color = 'black', linetype = 'dashed', alpha = 0.5) +
  ylab("Marginal Effect") +
  xlab(element_blank()) +
  ylim(-0.012, 0.012) +
  theme_bw() +
  theme(
    axis.title.x = element_text(size=14),
    axis.title.y = element_text(size=14),
    axis.text.y =  element_text(size=12),
    axis.text.x =  element_text(size=12)
  )




#
# ---- 4.3 Create Interaction plot for dominance and election proximity (Figure 6) -----
#

logit_int1 <- glm(subst_policy ~ 
                  proximity_nextelection*dominance_bin+
                  proximity_btw+
                  partei,
                data = groupappeals_fullset_v1 , 
                family= "binomial")


plot_model(logit_int1, 
           type = "pred", 
           terms = c("proximity_nextelection", "dominance_bin"), 
           ci.lvl = 0.95,
           show.data = FALSE) +
  scale_color_manual(values = c("black", "grey"), 
                     name = NULL,  # Remove legend title
                     labels = c("Challenger parties", "Dominant parties")) +  # Set custom labels
  scale_fill_manual(values = c("black", "grey"), 
                    guide = "none") +  # Remove fill legend if not needed
  theme_minimal() +
  xlab("Proximity to next election") +
  ylab("Propensity to use policy-based group appeal") +
  labs(title = NULL) +
  theme(legend.position = "right",
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())+
  theme_bw() +
  theme(
    axis.title.x = element_text(size=10),
    axis.text.y =  element_text(size=10),
    axis.text.x =  element_text(size=10)
  )













#
#---- 5. Appendix------
#

#
# ------ 5.1 Governing times per party (Figure A1) ----
#

ggplot(dominance_since1990,
       aes(y = land,
           yend = land,
           x = dominance_total/365,
           xend = 0)) +
  geom_segment() +
  geom_point(colour = "black") +
  theme_bw() +
  theme(legend.position = "none") +
  facet_wrap(~partei, scales = "free_y") +  
  labs(y = "", x = "") +
  theme(
    text = element_text(family = "Arial", angle = 0) 
  )

#
# ---- 5.2 Performance metrics of the social groups dictionary (Table A1) ----
#

# sample 800 tweets in order to check validity
# validity_dictionary <- raw_tweets_uncoded[sample(nrow(raw_tweets_uncoded), 
#                                                  800), ] # raw_tweets_uncoded
#                                                          # not loaded here due
#                                                          # to twitter policy
# # validity dictionary is exported and hand-coded for validity
# 
# # read coded cases (800 validated tweets)
# val_dict <- read.xlsx("Data/validity_dictionary.xlsx")
# 
# # Calculate Accuracy
# accuracy <- (sum(val_dict$TP) + sum(val_dict$TN)) /
#   (sum(val_dict$TP) + sum(val_dict$FP) + sum(val_dict$TN) + sum(val_dict$FN))
# 
# # Calculate Precision
# precision <- sum(val_dict$TP) / (sum(val_dict$TP) + sum(val_dict$FP))
# 
# # Calculate Recall
# recall <- sum(val_dict$TP) / (sum(val_dict$TP) + sum(val_dict$FN))
# 
# # Calculate F1
# f1_score <- (2 * precision * recall) / (precision + recall)

#
# ----- 5.3  Regression models with alternative dominance specifications (Table A3) -----
#
logit_m2_alt <- glm(subst_policy ~ dominance_wgt_year +
                  partei,
                  #as.factor(land),
                data = groupappeals_fullset_v1 , 
                family= "binomial")

logit_m4_alt <- glm(subst_policy ~ dominance_wgt_year*MP_wgt_year + 
                   partei + 
                   gov,
                   #as.factor(land),
                 data = groupappeals_fullset_v1, 
                 family= "binomial")


# combine models in one table and export as html
models <- list(logit_m2_alt, logit_m4_alt)

# Compute odds ratios, confidence intervals, and p-values
or_list <- lapply(models, extract_or)
ci_list <- lapply(models, compute_ci)
p_values_list <- lapply(models, extract_p_values)

# Generate the regression table using stargazer
stargazer(models, 
          type = "html",
          coef = or_list,  # Display odds ratios instead of log-odds
          ci = TRUE,  # Enable confidence intervals
          ci.custom = ci_list,  # Provide exponentiated confidence intervals
          p = p_values_list,  # Ensure correct p-values are used
          omit = "land|partei",  # Omit specified variables
          digits = 3,
          title = "Regression Table with Odds Ratios and Confidence Intervals",
          model.names = FALSE,  # No default model names
          column.labels = c("Model 2", "Model 6"),  # Custom model names
          dep.var.labels.include = FALSE,
          add.lines = list(
            c("State dummies", "YES", "YES"),
            c("Party dummies", "YES", "YES")
          ),
          out = "Figures/regtable_annex_weighted_mar2025.html")



#
# ---- 5.4 Additional model specifications (Table A5) -----
#

logit_a1 <- glm(subst_policy ~ dominance_bin +
                  #partei+
                  as.factor(land),
                data = groupappeals_fullset_v1 , 
                family= "binomial")


logit_a2 <- glm(subst_policy ~ dominance_year +
                  #partei+
                  as.factor(land),
                data = groupappeals_fullset_v1 , 
                family= "binomial")

logit_a3 <- glm(subst_policy ~ proximity_nextelection +
                  proximity_btw+
                  partei,
                #as.factor(land),
                data = groupappeals_fullset_v1 , 
                family= "binomial")

logit_a4 <- glm(subst_policy ~ proximity_nextelection +
                  proximity_btw+
                  #partei,
                  as.factor(land),
                data = groupappeals_fullset_v1 , 
                family= "binomial")

logit_a5 <- glm(subst_policy ~ MP_year +
                  partei,
                #as.factor(land),
                data = groupappeals_fullset_v1 , 
                family= "binomial")

logit_a6 <- glm(subst_policy ~ MP_year +
                  #partei,
                  as.factor(land),
                data = groupappeals_fullset_v1 , 
                family= "binomial")

logit_a7 <- glm(subst_policy ~ dominance_year*MP_year + gov +
                  partei+
                  as.factor(land),
                data = groupappeals_fullset_v1 , 
                family= "binomial")



#Define models
models <- list(logit_a1, logit_a2, logit_a3, logit_a4, logit_a5, logit_a6, logit_a7)

# Compute odds ratios, confidence intervals, and p-values
or_list <- lapply(models, extract_or)
ci_list <- lapply(models, compute_ci)
p_values_list <- lapply(models, extract_p_values)

# Generate the regression table using stargazer
stargazer(models, 
          type = "html",
          coef = or_list,  # Display odds ratios instead of log-odds
          ci = TRUE,  # Enable confidence intervals
          ci.custom = ci_list,  # Provide exponentiated confidence intervals
          p = p_values_list,  # Ensure correct p-values are used
          omit = "land|partei",  # Omit specified variables
          digits = 3,
          title = "Regression Table (Extra Appendix) with Odds Ratios and CIs",
          dep.var.labels.include = FALSE,
          add.lines = list(
            c("State dummies", "YES", "YES", "YES", "YES", "YES", "YES"),
            c("Party dummies", "NO", "NO", "YES", "YES", "YES", "YES")
          ),
          out = "Figures/regtable_extraappendix_mar2025.html")


#
# ---- 5.5  Policy-based appeals aggregated to national parties(Figure A2) -----
#

# create grouped df
byparty <- group_by(subst_perparty, partei) %>%
  summarise(
    count = n(),
    mean = mean(prop_subst, na.rm = TRUE),
    sd = sd(prop_subst, na.rm = TRUE)
  )

# change party variable to factor
byparty$partei <- as.factor(byparty$partei)

# create boxplot
ggboxplot(subst_perparty, x = "partei", y = "prop_subst",
                             ylab = "Proportion of substantial group appeals", xlab = "Party") +
  scale_y_continuous(limits = c(0,0.7), 
                     labels = c("0", "20", "40", "60"), 
                     breaks = c(0,0.2,0.4,0.6), 
                     name = "Share of policy-based group appeals") +
  theme_classic() +
  theme(
    axis.title.x = element_text(size = 14),
    axis.title.y = element_text(size = 14),
    axis.text.y =  element_text(size = 12),
    axis.text.x =  element_text(size = 12)
  )




#
# ---- 5.6 Regression models without all AfD branches (Table A7) -----
#
withoutafd <- groupappeals_fullset_v1 %>% filter(partei != "AfD")

# Calculate logit models without AfD
logit_woafd_m3 <- glm(subst_policy ~ dominance_bin +
                        partei,
                  #as.factor(land),
                data = withoutafd , 
                family = "binomial")

logit_woafd_m5 <- glm(subst_policy ~ dominance_year +
                  partei,
                  #as.factor(land),
                data = withoutafd , 
                family = "binomial")


logit_woafd_m7 <- glm(subst_policy ~ dominance_bin*proximity_nextelection+
                        proximity_btw+
                        partei,
                        #as.factor(land),
                      data = withoutafd , 
                      family = "binomial")

logit_woafd_m11 <- glm(subst_policy ~ dominance_year*MP_year +  partei+ gov,
                   #as.factor(land),
                 data = withoutafd , 
                 family = "binomial")


# combine models in one table with texreg

# Define models
models <- list(logit_woafd_m3, logit_woafd_m5, logit_woafd_m7, logit_woafd_m11)

# Compute odds ratios, confidence intervals, and p-values
or_list <- lapply(models, extract_or)
ci_list <- lapply(models, compute_ci)
p_values_list <- lapply(models, extract_p_values)

# Generate the regression table using stargazer
stargazer(models, 
          type = "html",
          coef = or_list,  # Display odds ratios instead of log-odds
          ci = TRUE,  # Enable confidence intervals
          ci.custom = ci_list,  # Provide exponentiated confidence intervals
          p = p_values_list,  # Ensure correct p-values are used
          omit = "land|partei",  # Omit specified variables
          digits = 3,
          title = "Regression Table (Excluding AfD) with Odds Ratios and CIs",
          dep.var.labels.include = FALSE,
          add.lines = list(
            c("Party dummies", "YES", "YES", "YES", "YES")
          ),
          out = "Figures/regtablex_wo_afd_mar2025.html")





# ---- 5.7 Dominance types (Figure A4)-----

# used to govern, never governed, permanent dominance, current government

dominance_since1990 <- dominance_since1990 %>%
  mutate(type = case_when(dominance_total == 0 ~ "never_governed",
                          before_2010 >= 0 & after_2010 == 0 & before_2000 > 0 ~ "usedtogovern",
                          before_2010 > 0 & after_2010 > 0 & before_2000 > 0 & MP_total > 0 ~ "permanently_dominant",
                          #dominance_total > 8000 ~ "permanently_dominant",
                          before_2010 == 0  & before_2000 == 0 & after_2010 > 0 ~ "newly_governing"))

# merge share of policy-based appeals to dominance data 
subst_perparty <- merge(subst_perparty, dominance_since1990, by = c("partei", 
                                                                    "land"))

# Filter out rows where type is NA
filtered_data <- subset(subst_perparty, !is.na(type))

#Calculate the median of prop_subst for each type
median_data <- subst_perparty %>%
  group_by(type) %>%
  summarise(median_prop_subst = median(prop_subst, na.rm = TRUE)) %>%
  arrange(median_prop_subst)  # Arrange by median prop_subst

# Reorder the levels of type based on median prop_subst
filtered_data$type <- factor(filtered_data$type, levels = median_data$type)

# Create a boxplot with reordered levels
ggplot(filtered_data, aes(x = type, y = prop_subst, fill = type)) +
  geom_boxplot() +
  stat_summary(fun=mean, geom="point", shape=17, size=3, color="black", position=position_dodge(0.75)) +  # Add means as triangular dots
  labs(title = "Proportion of policy-based group appeals per government type",
       x = "Party Type",
       y = "Proportion of policy-based group appeals") +
  theme_minimal() +
  theme(legend.position = "none")



#
# ---- 5.8 Tabulated data (Table A9) ----
#

# Create the table
table_data <- table(groupappeals_fullset_v1$dominance_bin, groupappeals_fullset_v1$subst_policy)

# Convert to percentages by row
percentages <- prop.table(table_data, margin = 1) * 100

# Print the result
percentages




